# File name: 	biv_chao_1998.R
# By: 		tw
# Purpose:		bayesian IV according to Chao, Philips (1998 JoE) with AdMit package
# Last Version:	07.02.09

setwd("D:/Nauka/_R")
library(foreign)
library(fBasics)
library(AdMit)

d <- read.dta("truth.dta")

# ivregress 2sls Y_IV (D_heli = Z_heli), first
# ivregress 2sls Y_IV_weak (D_heli_weak = Z_heli), first
# ivreg Y_IV (D_heli = Z_heli), first
# ivreg Y_IV_weak (D_heli_weak = Z_heli), first
# reg Y_IV_full Z_heli



'KERNEL' <- function(par,Y,D,Z,log=TRUE){
if (is.vector(par)) par <- t(par)

    'KERNEL.SUB' <- function(par,Y,D,Z){

	# logLikelihood function first
		v <- matrix(0,length(Y),2)
		v[,1] <- Y - par[3]*par[1] - par[2] - par[4]*par[1]*Z
		v[,2] <- D - par[3] - par[4]*Z

		Om <- matrix(0,2,2)
		Om[1,1] <- par[5]
		Om[2,1] <- par[6]
		Om[1,2] <- par[6]				
		Om[2,2] <- par[7]

		ll <- -(length(Y)/2)*log(abs(det(Om))) -0.5*tr(solve(Om)%*%t(v)%*%v)
		
	# log of Prior now
		Q <- matrix(-(1/length(Y)),length(Y),length(Y))
		for (i in 1:length(Y)) Q[1,1] <- Q[1,1]+1

		lp <- -1.5*log(abs(det(Om))) + .5*log(abs(det((par[4]^2)*t(Z)%*%Q%*%Z)))
	
		return(ll+lp)
    }
    r <- apply(par,1,KERNEL.SUB,Y=Y,D=D,Z=Z)
    if (!log) r <- exp(r)
    return(r)
}

p <- matrix(0,10000,7)
p[,1] <- runif(10000,0,.1)
p[,2] <- runif(10000,15,15.5)
p[,3] <- runif(10000,0,.2)
p[,4] <- runif(10000,.5,1)
p[,5] <- runif(10000,.5,3)
p[,6] <- runif(10000,0,.25)
p[,7] <- runif(10000,.5,3)

k <- KERNEL(p,d$Y_IV,d$D_heli,d$Z_heli)
print(k[which.max(k)])
par <- p[which.max(k),]
par <- c(0.02236423, 15.09986261, 0.09221081, 0.84413780, 1.13284599, 0.11078167, 0.50015638)
par <- c(.047306,15.17836,.102459,.7998847,1,.1,1)
par <- c(0.08166995, 15.08989957,  0.10522812,  0.69945933,  1.34873008,  0.14732127,0.50143253)

KERNEL(par,d$Y_IV,d$D_heli,d$Z_heli)

iv.admit <- AdMit(KERNEL,mu0=par,control=list(trace=TRUE,trace.mu=TRUE),Y=d$Y_IV,D=d$D_heli,Z=d$Z_heli)

'G' <- function(par){
	return(par)
}
iv.admitis <- AdMitIS(N=10000,KERNEL,G=G,mit=iv.admit$mit,
)

iv.admitmh <- AdMitMH(N=10000,KERNEL,mit=iv.admit$mit,Y=d$Y_IV,D=d$D_heli,Z=d$Z_heli)

plot(density(iv.admitmh$draws[,1]))
x11()
plot(density(iv.admitmh$draws[,2]))

save(mu0,iv.admit,iv.admitis,iv.admitmh,file="biv_lennart_2007_results.RData")